FigureYa292HCCsubtype

Author(s): Xiaofan Lu; Yasi Zhang
Reviewer(s): Ying Ge
Date: 2025-05-20

Academic Citation

If you use this code in your work or research, we kindly request that you cite our publication:

Xiaofan Lu, et al. (2025). FigureYa: A Standardized Visualization Framework for Enhancing Biomedical Data Interpretation and Research Efficiency. iMetaMed. https://doi.org/10.1002/imm3.70005

需求描述

Requirements Description

主要是红框中的图，都是公共数据和已发表的文章。

我们随后比较了基于TMDU和TCGA研究多平台数据的分类。

根据以往的分子分型marker预测现有样本属于哪一类，再跟当前分型比较，计算p值。

Mainly the figures in the red box, which are all public data and published articles.

We then compared our classification based on multi-platform data from the TMDU and TCGA studies.

Using previous molecular subtype markers to predict which category the existing samples belong to, and then comparing them with the current classification to calculate the p-value.

出自：https://linkinghub.elsevier.com/retrieve/pii/S2352396418306340

图3. 肝细胞癌分子分型概述。(a) TMDU测试研究（上）和TCGA验证研究（下）中，基于既往定义的肝细胞癌分子分型相关基因集的综合评分比较。(b) 分子亚型的示意图。

Source: https://linkinghub.elsevier.com/retrieve/pii/S2352396418306340

Fig. 3. Summary ofmolecular classification of HCC. (a) Comparison of aggregate scores with gene sets associated with the previously definedmolecular classifications of HCC in the TMDU test study (upper) and TCGA validation study (lower). (b) Schematic representation of molecular subtypes.

2017年的这篇Cell文章里就出现了类似的图：

A similar figure appeared in this 2017 Cell paper:

出自：https://linkinghub.elsevier.com/retrieve/pii/S0092867417306396

图2. 肝癌呈现独特的基因高甲基化模式。 (A) 基于肝癌组织相对于正常组织的基因高甲基化无监督聚类分析揭示了四个明显不同的亚群。研究人员分析了约15,000个在196名HCC患者中显示显著高甲基化的CpG位点，并以热图形式展示，其中正常组织和肿瘤样本按聚类分组排列在列中。每个CpG位点的甲基化强度通过行显示。热图上方标注了四个不同的高甲基化聚类群，下方条形图展示了各聚类群中个体肿瘤的临床和分子特征分布。右侧p值表示各特征的非随机分布具有统计学显著性。

Source: https://linkinghub.elsevier.com/retrieve/pii/S0092867417306396

Figure 2. Liver Cancers Show Distinct Gene Hypermethylation Patterns. (A) Unsupervised clustering analysis of gene hypermethylation in HCC tumor relative to normal tissue reveals four distinct subgroups. Roughly 15,000 CpG sites showing significant hypermethylation in 196 HCC patients were analyzed and are shown in heatmap format with normal tissues and tumors organized in columns according to cluster designation. Intensity of methylation for each CpG site is indicated by row. Above the heatmap the four distinct hypermethylation clusters are shown, and below are bars indicating the distribution of clinical and molecular attributes of the individual tumors by cluster. To the right, p values indicate significant non-random distributions for each attribute.

应用场景

Application Scenario

根据表达谱对HCC样本进行聚类，并和其他已有的分类器做比较。

Clustering of HCC samples based on expression profiles and comparison with other existing classifiers.

环境设置

Environment Setup

source("install_dependencies.R")

## Starting R package installation...
## ===========================================
## 在 Linux 系统上，建议安装以下系统依赖:
## sudo apt-get install -y libcurl4-openssl-dev libssl-dev libxml2-dev
## 
## Installing CRAN packages...
## Package already installed: ClassDiscovery 
## Package already installed: gplots 
## Package already installed: tidyverse 
## 
## Installing Bioconductor packages...
## Package already installed: CMScaller 
## Package already installed: ComplexHeatmap 
## Package already installed: clusterProfiler 
## Package already installed: org.Hs.eg.db 
## 
## ===========================================
## Package installation completed!
## 
## 验证包安装状态:
## ✓ ClassDiscovery is installed
## ✓ gplots is installed
## ✓ tidyverse is installed
## ✓ CMScaller is installed
## ✓ ComplexHeatmap is installed
## ✓ clusterProfiler is installed
## ✓ org.Hs.eg.db is installed
## You can now run your R scripts in this directory.

# 用来实现nearest template prediction(NTP)
# Implementation of Nearest Template Prediction (NTP)
library(CMScaller) 

## CMScaller v2.0.1; genome annotation: GENCODE v32/GRCh38.p13

library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.2     ✔ tibble    3.3.0
## ✔ lubridate 1.9.4     ✔ tidyr     1.3.1
## ✔ purrr     1.0.4

## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(clusterProfiler)

## 
## clusterProfiler v4.14.6 Learn more at https://yulab-smu.top/contribution-knowledge-mining/
## 
## Please cite:
## 
## T Wu, E Hu, S Xu, M Chen, P Guo, Z Dai, T Feng, L Zhou, W Tang, L Zhan,
## X Fu, S Liu, X Bo, and G Yu. clusterProfiler 4.0: A universal
## enrichment tool for interpreting omics data. The Innovation. 2021,
## 2(3):100141
## 
## Attaching package: 'clusterProfiler'
## 
## The following object is masked from 'package:purrr':
## 
##     simplify
## 
## The following object is masked from 'package:stats':
## 
##     filter

library(org.Hs.eg.db)

## Loading required package: AnnotationDbi
## Loading required package: stats4
## Loading required package: BiocGenerics
## 
## Attaching package: 'BiocGenerics'
## 
## The following objects are masked from 'package:lubridate':
## 
##     intersect, setdiff, union
## 
## The following objects are masked from 'package:dplyr':
## 
##     combine, intersect, setdiff, union
## 
## The following objects are masked from 'package:stats':
## 
##     IQR, mad, sd, var, xtabs
## 
## The following objects are masked from 'package:base':
## 
##     anyDuplicated, aperm, append, as.data.frame, basename, cbind,
##     colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,
##     get, grep, grepl, intersect, is.unsorted, lapply, Map, mapply,
##     match, mget, order, paste, pmax, pmax.int, pmin, pmin.int,
##     Position, rank, rbind, Reduce, rownames, sapply, saveRDS, setdiff,
##     table, tapply, union, unique, unsplit, which.max, which.min
## 
## Loading required package: Biobase
## Welcome to Bioconductor
## 
##     Vignettes contain introductory material; view with
##     'browseVignettes()'. To cite Bioconductor, see
##     'citation("Biobase")', and for packages 'citation("pkgname")'.
## 
## Loading required package: IRanges
## Loading required package: S4Vectors
## 
## Attaching package: 'S4Vectors'
## 
## The following object is masked from 'package:clusterProfiler':
## 
##     rename
## 
## The following objects are masked from 'package:lubridate':
## 
##     second, second<-
## 
## The following objects are masked from 'package:dplyr':
## 
##     first, rename
## 
## The following object is masked from 'package:tidyr':
## 
##     expand
## 
## The following object is masked from 'package:utils':
## 
##     findMatches
## 
## The following objects are masked from 'package:base':
## 
##     expand.grid, I, unname
## 
## 
## Attaching package: 'IRanges'
## 
## The following object is masked from 'package:clusterProfiler':
## 
##     slice
## 
## The following object is masked from 'package:lubridate':
## 
##     %within%
## 
## The following objects are masked from 'package:dplyr':
## 
##     collapse, desc, slice
## 
## The following object is masked from 'package:purrr':
## 
##     reduce
## 
## 
## Attaching package: 'AnnotationDbi'
## 
## The following object is masked from 'package:clusterProfiler':
## 
##     select
## 
## The following object is masked from 'package:dplyr':
## 
##     select

library(ClassDiscovery)

## Loading required package: cluster
## Loading required package: oompaBase

library(ComplexHeatmap)

## Loading required package: grid
## ========================================
## ComplexHeatmap version 2.25.1
## Bioconductor page: http://bioconductor.org/packages/ComplexHeatmap/
## Github page: https://github.com/jokergoo/ComplexHeatmap
## Documentation: http://jokergoo.github.io/ComplexHeatmap-reference
## 
## If you use it in published research, please cite either one:
## - Gu, Z. Complex Heatmap Visualization. iMeta 2022.
## - Gu, Z. Complex heatmaps reveal patterns and correlations in multidimensional 
##     genomic data. Bioinformatics 2016.
## 
## 
## The new InteractiveComplexHeatmap package can directly export static 
## complex heatmaps into an interactive Shiny app with zero effort. Have a try!
## 
## This message can be suppressed by:
##   suppressPackageStartupMessages(library(ComplexHeatmap))
## ========================================

library(gplots)

## 
## Attaching package: 'gplots'
## 
## The following object is masked from 'package:oompaBase':
## 
##     redgreen
## 
## The following object is masked from 'package:IRanges':
## 
##     space
## 
## The following object is masked from 'package:S4Vectors':
## 
##     space
## 
## The following object is masked from 'package:stats':
## 
##     lowess

# 显示英文报错信息
# Show English error messages
Sys.setenv(LANGUAGE = "en") 

# 禁止chr转成factor
# Prevent character to factor conversion
options(stringsAsFactors = FALSE) 

自定义函数

Custom Functions

standarize.fun <- function(indata=NULL, halfwidth=NULL, centerFlag=T, scaleFlag=T) {  
  outdata=t(scale(t(indata), center=centerFlag, scale=scaleFlag))
  if (!is.null(halfwidth)) {
    outdata[outdata>halfwidth]=halfwidth
    outdata[outdata<(-halfwidth)]= -halfwidth
  }
  return(outdata)
}

fpkmToTpm <- function(fpkm)
{
  exp(log(fpkm) - log(sum(fpkm)) + log(1e6))
}

输入文件

Input Files

data_mutations.txt，突变数据，下载自cBioPortalhttps://www.cbioportal.org/。

LIHC.htseq_fpkm.tsv.gz，表达谱数据FPKM，已经过log2(fpkm+1)转换，下载地址：https://xenabrowser.net/datapages/?dataset=TCGA-LIHC.htseq_fpkm.tsv&host=https%3A%2F%2Fgdc.xenahubs.net&removeHub=https%3A%2F%2Fxena.treehouse.gi.ucsc.edu%3A443

gencode.v22.annotation.gene.probeMap，ID/Gene Mapping，下载地址同上。

data_mutations.txt，mutation data, downloaded from cBioPortalhttps://www.cbioportal.org/.

TCGA-LIHC.star_fpkm.tsv.gz, expression profile data FPKM, already converted by log2(fpkm+1), download link: https://xenabrowser.net/datapages/?dataset=TCGA-LIHC.star_fpkm.tsv&host=https%3A%2F%2Fgdc.xenahubs.net&removeHub=https%3A%2F%2Fxena.treehouse.gi.ucsc.edu%3A443

gencode.v36.annotation.gtf.gene.probemap, ID/Gene Mapping, download link same as above.

# 读取突变数据
# Read mutation data
maf <- read_tsv("data_mutations.txt", comment = "#")

## Rows: 223 Columns: 114
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (100): Hugo_Symbol, Center, NCBI_Build, Strand, Consequence, Variant_Cla...
## dbl  (14): Entrez_Gene_Id, Chromosome, Start_Position, End_Position, t_ref_c...
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

# 把突变数据转成01矩阵（方法跟FigureYa288MutualExclusivity一样）
# Convert mutation data to binary matrix (same method as FigureYa288MutualExclusivity)
mut.binary <- matrix(0,nrow = length(unique(maf$Hugo_Symbol)),ncol = length(unique(maf$Tumor_Sample_Barcode)),dimnames = list(unique(maf$Hugo_Symbol),unique(maf$Tumor_Sample_Barcode)))

# 填充突变矩阵(1表示存在突变)
# Fill in mutation matrix (1=mutation present) 
for (i in colnames(mut.binary)) {
  tmp <- maf[which(maf$Tumor_Sample_Barcode == i),]
  tmp <- tmp[which(tmp$Variant_Classification %in% c("Frame_Shift_Del", "Frame_Shift_Ins", "Splice_Site", "Translation_Start_Site","Nonsense_Mutation", "Nonstop_Mutation", "In_Frame_Del","In_Frame_Ins", "Missense_Mutation")),]
  for (j in tmp$Hugo_Symbol)
    mut.binary[j,i] <- 1
}

# 转为数据框
# Convert to dataframe 
mut.binary <- as.data.frame(mut.binary)

# 读取表达谱数据
# Read expression profiles
fpkm <- read.delim("TCGA-LIHC.star_fpkm.tsv.gz", sep = "\t", row.names = 1, check.names = F, stringsAsFactors = F, header = T)
fpkm <- 2^fpkm - 1

# 把fpkms转为tpm
# Convert FPKM to TPM
tpm <- apply(fpkm, 2, fpkmToTpm)
tpm <- as.data.frame(log2(tpm + 1))
rm(fpkm); gc()

##            used  (Mb) gc trigger   (Mb)  max used   (Mb)
## Ncells  6188058 330.5   11557434  617.3   8171455  436.5
## Vcells 37404520 285.4  193310788 1474.9 241605833 1843.4

# 加载基因注释
# Load gene annotation
Ginfo <- read.delim("gencode.v36.annotation.gtf.gene.probemap",row.names = 1,sep = "\t",check.names = F,stringsAsFactors = F,header = T)

# 匹配表达谱和注释的基因
# Match genes between expression and annotation
comgene <- intersect(rownames(Ginfo), rownames(tpm))
Ginfo <- Ginfo[comgene,]
tpm <- tpm[comgene,]
identical(rownames(tpm), rownames(Ginfo))

## [1] TRUE

# 按基因名聚合取表达中位数
# Aggregate duplicate genes by median expression
tpm$Gene <- Ginfo[rownames(tpm),"gene"]
tpm <- as.data.frame(apply(tpm[,setdiff(colnames(tpm), "Gene")], 2, function(x) tapply(x, INDEX = factor(tpm$Gene), FUN=median, na.rm = TRUE))) 

# 截短样本ID
# Trim sample IDs
colnames(tpm) <- substr(colnames(tpm),1,15)

expr <- tpm

# 提取突变和表达的共同样本
# Extract common samples between mutation and expression
comsam <- intersect(colnames(expr), colnames(mut.binary))
mut.binary <- mut.binary[,comsam]
expr <- expr[,comsam]

# 加载原文用于聚类的差异表达基因
# Load DEGs for clustering from original paper
degs <- read.table("degs.txt",sep = "\t",row.names = NULL,header = T,check.names = F,stringsAsFactors = F)
comgene <- intersect(degs$Gene, rownames(expr))

# 从论文中获取iCluster亚型：肝细胞癌的综合基因组特征分析
# Get iCluster subtypes from paper: Comprehensive and Integrative Genomic Characterization of Hepatocellular Carcinoma
icluster <- read.delim("icluster.txt",sep = "\t",row.names = 2,check.names = F,stringsAsFactors = F,header = T)
rownames(icluster) <- substr(rownames(icluster),1,15)

Nearest template prediction

加载不同HCC亚型签名并制作nearest template prediction所需的模版。

Loading HCC subtype signatures and creating templates for nearest template prediction.

boyault <- read.delim("Boyault-marker-all.txt",sep = "\t",row.names = NULL,check.names = F,stringsAsFactors = F,header = T)

# 将Entrez ID转换为基因符号
# Convert Entrez IDs to gene symbols
tmp <- bitr(boyault$`Gene ID #1`,fromType = "ENTREZID",toType = "SYMBOL",OrgDb = org.Hs.eg.db)

## 'select()' returned 1:1 mapping between keys and columns

## Warning in bitr(boyault$`Gene ID #1`, fromType = "ENTREZID", toType = "SYMBOL",
## : 0.31% of input gene IDs are fail to map...

# 提取关键列
# Extract relevant columns 
boyault <- boyault[,c(2,4)]
colnames(boyault) <- c("ENTREZID","class")

# 合并基因符号映射
# Merge with symbol mapping 
boyault <- merge(tmp,boyault,by = "ENTREZID", all.x = T)
colnames(boyault)[2] <- "probe"

chiang <- read.delim("Chiang-marker-all.txt",sep = "\t",row.names = NULL,check.names = F,stringsAsFactors = F,header = T)
chiang <- chiang[,c(1,4)]
colnames(chiang) <- c("probe","class")

hoshida <- read.delim("Hoshida-genelist.txt",sep = "\t",row.names = NULL,check.names = F,stringsAsFactors = F,header = T)
hoshida <- hoshida[,c(1,2)]
colnames(hoshida) <- c("probe","class")
hoshida$class <- paste0("C",hoshida$class)

lee <- read.delim("Lee-marker-all.txt",sep = "\t",row.names = NULL,check.names = F,stringsAsFactors = F,header = T)
lee <- lee[,c(1,3)]
colnames(lee) <- c("probe","class")

通过nearest template prediction获取其他亚型。

Obtaining other subtypes via nearest template prediction.

tcga.boyault <- ntp(emat = t(scale(t(expr))),
                    templates = boyault,
                    doPlot = T,
                    seed = 19991018)

tcga.chiang <- ntp(emat = t(scale(t(expr))),
                   templates = chiang,
                   doPlot = T,
                   seed = 19991018)

tcga.lee <- ntp(emat = t(scale(t(expr))),
                   templates = lee,
                   doPlot = T,
                   seed = 19991018)

tcga.hoshida <- ntp(emat = t(scale(t(expr))),
                    templates = hoshida,
                    doPlot = T,
                    seed = 19991018)

开始画图

Plotting

# 构建样本注释
# Construct sample annotation
annCol <- data.frame(CTNNB1 = ifelse(as.numeric(mut.binary["CTNNB1",]) == 0,"WT","MT"),
                     Lee = as.character(tcga.lee$prediction),
                     Hoshida = as.character(tcga.hoshida$prediction),
                     Boyault = as.character(tcga.boyault$prediction),
                     Chiang = as.character(tcga.chiang$prediction),
                     row.names = colnames(mut.binary),
                     stringsAsFactors = F)

# 整合iCluster亚型
# Integrate iCluster subtypes
annCol[intersect(rownames(annCol),rownames(icluster)),"iCluster"] <- icluster[intersect(rownames(annCol),rownames(icluster)), "iCluster clusters (k=3, Ronglai Shen)"]

# 标记缺失数据
# Mark missing data
annCol[is.na(annCol$iCluster),"iCluster"] <- "N/A"

# 标准化亚型标签
# Standardize subtype labels
annCol$Hoshida <- gsub("C","S",annCol$Hoshida)
annCol$iCluster <- gsub("iCluster:","iC",annCol$iCluster)

# 优化Chiang亚型标签
# Refine Chiang subtype labels
annCol[which(annCol$Chiang == "interferon class"), "Chiang"] <- "INTERFERON"
annCol[which(annCol$Chiang == "proliferation class"), "Chiang"] <- "PROLIFERATION"
annCol[which(annCol$Chiang == "CTNNB1 class"), "Chiang"] <- "CTNNB1"
annCol[which(annCol$Chiang == "unannotated class"), "Chiang"] <- "UNANNOTATED"
annCol[which(annCol$Chiang == "chromosome 7 polysomy class"), "Chiang"] <- "POLYSOMY7"

# 设置注释列的因子水平
# Set factor levels for annotation columns
annCol$CTNNB1 <- factor(annCol$CTNNB1, levels = c("MT","WT"))
annCol$Lee <- factor(annCol$Lee, levels = c("SURVIVAL_DN","SURVIVAL_UP"))
annCol$Hoshida <- factor(annCol$Hoshida, levels = c("S1","S2","S3"))
annCol$Boyault <- factor(annCol$Boyault, levels = c("G12","G3","G56"))
annCol$Chiang <- factor(annCol$Chiang, levels = c("PROLIFERATION","CTNNB1","INTERFERON","POLYSOMY7","UNANNOTATED"))
annCol$iCluster <- factor(annCol$iCluster, levels = c("iC1","iC2","iC3","N/A"))

# 创建注释颜色方案
# Create color scheme for annotations
annColors <- list()
annColors[["CTNNB1"]] <- c("WT" = "white","MT" = "black")
annColors[["Lee"]] <- c("SURVIVAL_DN" = "#E80035","SURVIVAL_UP" = "#DBDEDD")
annColors[["Hoshida"]] <- c("S1" = "#E8536B","S2" = "#F6B879","S3" = "#DFEBAF")
annColors[["Boyault"]] <- c("G12" = "#E8536B","G3" = "#F6B879","G56" = "#DFEBAF")
annColors[["Chiang"]] <- c("PROLIFERATION" = "#F2A1A2","CTNNB1" = "#68BE8B","INTERFERON" = "#2CA8E1","POLYSOMY7" = "#A2D9F1","UNANNOTATED" = "#BAC8E5")
annColors[["iCluster"]] <- c("iC1" = "#A8BFDE","iC2" = "#DA8A88","iC3" = "#BCA3D6","N/A" = "white")
annColors[["Group"]] <- c("A" = "#FE0000","B" = "black")
annColors[["MS"]] <- c("MS1" = "#FE0000", "MS2" = "#00AF50","MS3" = "#0071C0")

# 进行无监督聚类
# Perform unsupervised clustering
indata <- expr[comgene,]
indata <- indata[rowSums(indata) > 0,]

# 计算距离矩阵
# Calculate distance matrices
hcs <- hclust(distanceMatrix(as.matrix(indata), "euclidean"), "ward.D")
hcg <- hclust(distanceMatrix(as.matrix(t(indata)), "euclidean"), "ward.D")

# 将树切割为3类
# Cut tree into 3 clusters
group <- cutree(hcs, k = 3)

# 检查样本分布和CTNNB1突变情况
# Check sample distribution and CTNNB1 mutation
table(group, annCol$CTNNB1) 

##      
## group MT WT
##     1 48  0
##     2 21 37
##     3 25 52

# 分配亚型标签
# Assign subtype labels
annCol$MS <- ifelse(group == 1, "MS2", ifelse(group == 3, "MS1","MS3")) # 由于1组有更多CTNNB1，所以是MS2，由于3组样本更少，所以是MS1
annCol$Group <- ifelse(group == 3, "A", "B") # 由于3组是MS1所以为Group A

# 准备热图数据
# Prepare data for heatmap
plotdata <- standarize.fun(indata, halfwidth = 2) # 数据标准化用于绘图

# 生成热图
# Generate heatmap
hm <- pheatmap(plotdata,
         color = bluered(64),
         border_color = NA,
         cluster_rows = hcg,
         cluster_cols = hcs,
         treeheight_row = 30,
         treeheight_col = 30,
         #cutree_cols = 3,
         show_rownames = F,
         show_colnames = F,
         annotation_col = annCol[,c("iCluster","Chiang","Boyault","Hoshida","Lee","CTNNB1","MS","Group")],
         annotation_colors = annColors,
         cellwidth = 0.6,
         cellheight = 0.2)

# 保存热图
# Save heatmap
pdf("heatmap.pdf", width = 8,height = 6)
draw(hm, heatmap_legend_side = "left", annotation_legend_side = "left")
invisible(dev.off())

计算group和其他类型的独立性

Testing independence between groups and other types

# 初始化输出表格
# Initialize output tables
outTabGroup <- outTabMS <- NULL

# 检验Group/MS与其他分类的独立性
# Perform independence tests between Group/MS and other classifications
for (i in c("iCluster","Chiang","Boyault","Hoshida","Lee","CTNNB1")) {
  tmp <- annCol[,c(i,"Group")]
  set.seed(123)
  p <- fisher.test(table(tmp[,1],tmp[,2]), 
                   simulate.p.value = T)$p.value 
  
  outTabGroup <- rbind.data.frame(outTabGroup,
                                  data.frame(VarA = "Group",
                                             VarB = i,
                                             p = p),
                                  stringsAsFactors = F)
  
  # 检验MS与当前分类的关系
  # Test MS vs current classification 
  tmp <- annCol[,c(i,"MS")]
  set.seed(123)
  p <- fisher.test(table(tmp[,1],tmp[,2]),simulate.p.value = T)$p.value
  outTabMS <- rbind.data.frame(outTabMS,
                               data.frame(VarA = "MS",
                                          VarB = i,
                                          p = p),
                               stringsAsFactors = F)
}

# 把p-values保存到文件
# Save p-values to files
write.table(outTabGroup, file = "independent test between group and other classification.txt",sep = "\t",row.names = F,col.names = T,quote = F)
write.table(outTabMS, file = "independent test between ms and other classification.txt",sep = "\t",row.names = F,col.names = T,quote = F)

# 保存工作镜像
# Save workspace image
#save.image(file = "LIHC.RData")

后期加工

Post-processing

输出的PDF文件是矢量图，可以用Illustrator等矢量图编辑工具打开，添加p value。如果想用代码添加p-values，可参考FigureYa165heatmapPvalue或FigureYa280TMEofSTS

The output PDF files are vector graphics that can be opened and edited with vector graphic tools like Illustrator to add p-values. If you prefer to add p-values programmatically, please refer to FigureYa165heatmapPvalue or FigureYa280TMEofSTS.

Session Info

sessionInfo()

## R version 4.4.1 (2024-06-14)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 20.04.4 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0 
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## time zone: Asia/Shanghai
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] grid      stats4    stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] gplots_3.2.0           ComplexHeatmap_2.25.1  ClassDiscovery_3.4.8  
##  [4] oompaBase_3.2.10       cluster_2.1.6          org.Hs.eg.db_3.20.0   
##  [7] AnnotationDbi_1.68.0   IRanges_2.40.1         S4Vectors_0.44.0      
## [10] Biobase_2.66.0         BiocGenerics_0.52.0    clusterProfiler_4.14.6
## [13] lubridate_1.9.4        forcats_1.0.0          stringr_1.5.1         
## [16] dplyr_1.1.4            purrr_1.0.4            readr_2.1.5           
## [19] tidyr_1.3.1            tibble_3.3.0           ggplot2_3.5.2         
## [22] tidyverse_2.0.0        CMScaller_2.0.1       
## 
## loaded via a namespace (and not attached):
##   [1] RColorBrewer_1.1-3      rstudioapi_0.17.1       jsonlite_2.0.0         
##   [4] shape_1.4.6.1           magrittr_2.0.3          magick_2.8.7           
##   [7] ggtangle_0.0.6          farver_2.1.2            rmarkdown_2.29         
##  [10] GlobalOptions_0.1.2     fs_1.6.6                zlibbioc_1.52.0        
##  [13] vctrs_0.6.5             Cairo_1.6-2             memoise_2.0.1          
##  [16] ggtree_3.14.0           htmltools_0.5.8.1       gridGraphics_0.5-1     
##  [19] sass_0.4.10             KernSmooth_2.23-24      bslib_0.9.0            
##  [22] plyr_1.8.9              cachem_1.1.0            igraph_2.0.3           
##  [25] lifecycle_1.0.4         iterators_1.0.14        pkgconfig_2.0.3        
##  [28] Matrix_1.7-3            R6_2.6.1                fastmap_1.2.0          
##  [31] gson_0.1.0              GenomeInfoDbData_1.2.13 clue_0.3-66            
##  [34] digest_0.6.37           aplot_0.2.5             enrichplot_1.26.6      
##  [37] colorspace_2.1-1        patchwork_1.3.0         RSQLite_2.4.1          
##  [40] timechange_0.3.0        httr_1.4.7              compiler_4.4.1         
##  [43] bit64_4.6.0-1           withr_3.0.2             doParallel_1.0.17      
##  [46] BiocParallel_1.40.2     DBI_1.2.3               R.utils_2.13.0         
##  [49] rjson_0.2.23            gtools_3.9.5            caTools_1.18.3         
##  [52] tools_4.4.1             ape_5.8-1               R.oo_1.27.1            
##  [55] glue_1.8.0              nlme_3.1-165            GOSemSim_2.32.0        
##  [58] reshape2_1.4.4          oompaData_3.1.5         fgsea_1.32.4           
##  [61] generics_0.1.4          gtable_0.3.6            tzdb_0.5.0             
##  [64] R.methodsS3_1.8.2       data.table_1.17.4       hms_1.1.3              
##  [67] XVector_0.46.0          ggrepel_0.9.6           foreach_1.5.2          
##  [70] pillar_1.10.2           yulab.utils_0.2.0       vroom_1.6.5            
##  [73] circlize_0.4.16         splines_4.4.1           treeio_1.30.0          
##  [76] lattice_0.22-5          bit_4.6.0               tidyselect_1.2.1       
##  [79] GO.db_3.20.0            Biostrings_2.74.1       knitr_1.50             
##  [82] xfun_0.52               matrixStats_1.5.0       stringi_1.8.7          
##  [85] UCSC.utils_1.2.0        lazyeval_0.2.2          ggfun_0.1.8            
##  [88] yaml_2.3.10             evaluate_1.0.3          codetools_0.2-19       
##  [91] qvalue_2.38.0           ggplotify_0.1.2         cli_3.6.5              
##  [94] jquerylib_0.1.4         dichromat_2.0-0.1       Rcpp_1.0.14            
##  [97] GenomeInfoDb_1.42.3     png_0.1-8               parallel_4.4.1         
## [100] blob_1.2.4              mclust_6.1.1            DOSE_4.0.1             
## [103] bitops_1.0-9            tidytree_0.4.6          scales_1.4.0           
## [106] crayon_1.5.3            GetoptLong_1.0.5        rlang_1.1.6            
## [109] cowplot_1.1.3           fastmatch_1.1-6         KEGGREST_1.46.0